#include "headers.hh"
#include "style.h"

const Int_t fColBl = kBlue;//7; //

void drawMeanGrSys01()
{
  TString tit;
  const TString pwd = gSystem->pwd();

  if(pwd.Contains("MC0")){
    tit+="ALICE Preliminary";
  }
  else{
    tit+="MC";
  }

  if(pwd.Contains("JPT5-10")){
    tit+=", #it{p}_{T,jet}^{ch} 5-10 GeV/#it{c}";
  }
  else if(pwd.Contains("JPT10-15")){
    tit+=", #it{p}_{T,jet}^{ch} 10-15 GeV/#it{c}";
  }
  else if(pwd.Contains("JPT15-20")){
    tit+=", #it{p}_{T,jet}^{ch} 15-20 GeV/#it{c}";
  }
  else if(pwd.Contains("JPT-5")){
    tit+=", inclusive";
  }
  else if(pwd.Contains("JPT20-100")){
    tit+=", #it{p}_{T,jet}^{ch} 20-100 GeV/#it{c}";
  }

  /*
  if(pwd.Contains("ETA0")){
    tit+=", | #eta^{track} | <0.2";
  }
  else if(pwd.Contains("ETA1")){
    tit+=", 0.2< | #eta^{track} | <0.4";
  }
  else if(pwd.Contains("ETA2")){
    tit+=", 0.4< | #eta^{track} | <0.6";
  }
  */
  //=======================================

  const Int_t col[]={fColBl, kRed, kGreen+3, kMagenta};//fColBl, fColBl, fColBl};
                     //

  TCanvas *c1 = new TCanvas;
  style::PadSetup(c1);
  c1->SetRightMargin(0.12);
  c1->SetTopMargin(0.1);
  gStyle->SetOptTitle(1);
  c1->SetBottomMargin(0.14);
  c1->SetLeftMargin(0.12);

  TFile::Open("outgrsys_id00.root");
  TGraphAsymmErrors *gFrac0=(TGraphAsymmErrors *)gDirectory->Get("grsys0");
  TGraphAsymmErrors *gFrac1=(TGraphAsymmErrors *)gDirectory->Get("grsys1");
  TGraphAsymmErrors *gFrac2=(TGraphAsymmErrors *)gDirectory->Get("grsys2");
  TGraphAsymmErrors *gFrac3=(TGraphAsymmErrors *)gDirectory->Get("grsys3");

  TFile::Open("outgrsys_id01.root");
  TGraphAsymmErrors *gMean0=(TGraphAsymmErrors *)gDirectory->Get("grsys0");
  TGraphAsymmErrors *gMean1=(TGraphAsymmErrors *)gDirectory->Get("grsys1");
  TGraphAsymmErrors *gMean2=(TGraphAsymmErrors *)gDirectory->Get("grsys2");
  TGraphAsymmErrors *gMean3=(TGraphAsymmErrors *)gDirectory->Get("grsys3");

  TFile::Open("outgrsys_id02.root");
  TGraphAsymmErrors *gReso0=(TGraphAsymmErrors *)gDirectory->Get("grsys0");
  TGraphAsymmErrors *gReso1=(TGraphAsymmErrors *)gDirectory->Get("grsys1");
  TGraphAsymmErrors *gReso2=(TGraphAsymmErrors *)gDirectory->Get("grsys2");
  TGraphAsymmErrors *gReso3=(TGraphAsymmErrors *)gDirectory->Get("grsys3");

  TGraphAsymmErrors * gFracs[]={gFrac0,gFrac1,gFrac2,gFrac3};
  TGraphAsymmErrors * gMeans[]={gMean0,gMean1,gMean2,gMean3};
  TGraphAsymmErrors * gResos[]={gReso0,gReso1,gReso2,gReso3};

  //============

  TFile *ff=new TFile("hist.root");
  TH2D *horiginal=(TH2D*)ff->Get("horiginal");
  horiginal->RebinX(5);
  horiginal->RebinY(10);

  TH2D * kk = (TH2D*)ff->Get("kk");
  TH1D *hyield = kk->ProjectionX();
  hyield->RebinX(5);

  gStyle->SetOptStat("");
  //  gSystem->Load("../../../../../../../tcf/lib/libstyle.so");
  style::ToNaturalScale(horiginal);
  horiginal->SetMinimum(0.9);

  style::ResetStyle(horiginal);

  horiginal->GetXaxis()->SetMoreLogLabels();

  horiginal->SetTitle(tit);
  gStyle->SetTitleY(0.975);
  gStyle->SetTitleW(0.95);
  horiginal->Draw("colz"); gPad->SetLogz();
  horiginal->SetXTitle("#it{p}^{track} (GeV/#it{c})");
  horiginal->SetYTitle("TPC signal (arb. unit)");
  horiginal->GetXaxis()->SetRangeUser(0.15,20);
  horiginal->GetYaxis()->SetRangeUser(20,150);
  
  TLatex * lco = new TLatex(3, 140, "TPC Coherent Fit");
  lco->Draw();

  TLatex *lpi=new TLatex(0.2, 55, "#pi"); //style::ResetStyle(lpi);
  lpi->Draw();

  TLatex *lp=new TLatex(0.55, 115, "p"); //style::ResetStyle(lp);
  lp->Draw();

  TLatex *lk=new TLatex(0.25, 120, "K"); //style::ResetStyle(lk);
  lk->Draw();

  TLatex *le=new TLatex(0.35, 76, "e"); //style::ResetStyle(le);
  le->Draw();

  for(int ii=0; ii<4; ii++){ 
    gMeans[ii]->SetFillColor(col[ii]); 
    gMeans[ii]->SetLineColor(col[ii]); 
    gMeans[ii]->SetFillStyle(0); 
    gMeans[ii]->SetLineWidth(2); 
    gMeans[ii]->SetLineColor(col[ii]); 
    gMeans[ii]->SetMarkerColor(col[ii]); 
    style::ResetStyle(gMeans[ii]);
    gMeans[ii]->Draw("pl");
  }
  gPad->SetLogx();
  gPad->SetLogz();

  c1->Print("data_mean.eps");
  c1->Print("data_mean.root");

  gPad->SetLogx(0);
  gPad->SetLogy(0);

  //-------

  for(Int_t ib=1; ib<= horiginal->GetNbinsX(); ib++){
    const Double_t cc = horiginal->Integral(ib, ib, 0, horiginal->GetNbinsY()+1);
    if(cc<50){
      continue;
    }

    const Double_t pp = pow(10, log10(horiginal->GetXaxis()->GetBinUpEdge(ib)*horiginal->GetXaxis()->GetBinLowEdge(ib))/2.);

    Double_t fitFrac[] = {-999, -999, -999, -999};
    Double_t fitErrFrac[] = {-999, -999, -999, -999};
    Double_t fitMean[] = {-999, -999, -999, -999};
    Double_t fitReso[] = {-999, -999, -999, -999};
    Double_t fitErrMean[] = {-999, -999, -999, -999};
    Double_t fitErrReso[] = {-999, -999, -999, -999};

    for(Int_t itype=0; itype<4; itype++){
      for(Int_t igr=0; igr<gMeans[itype]->GetN(); igr++){
        if(fabs(pp-gMeans[itype]->GetX()[igr])<EPSILON){
          fitMean[itype] = gMeans[itype]->GetY()[igr];
          fitErrMean[itype] = gMeans[itype]->GetEYlow()[igr];

          fitFrac[itype] = gFracs[itype]->GetY()[igr];
          fitErrFrac[itype] = gFracs[itype]->GetEYlow()[igr];

          fitReso[itype] = gResos[itype]->GetY()[igr];
          fitErrReso[itype] = gResos[itype]->GetEYlow()[igr];
        }
      }
    }
   
    if(fitMean[0]<-990){
      //pPb: ib 3 paritally exist  and the actual starting ib is 4
      for(Int_t igr=0; igr<gMeans[0]->GetN(); igr++){
        printf("continue grx not found! ib %d pp %e grx %e\n", ib, pp, gMeans[0]->GetX()[igr]);
      }
      continue;
    }


    TH1D * hb = horiginal->ProjectionY(Form("horiginal%d", ib), ib, ib);    
    //hb->SetTitle(Form("%s, #it{p}^{track} %.1f GeV/#it{c};(d#it{E}/d#it{x}-d#it{E}/d#it{x}_{#pi})/#sigma_{#pi};counts", tit.Data(), pp));
    hb->SetTitle(Form("%s;(d#it{E}/d#it{x}-d#it{E}/d#it{x}_{#pi})/#sigma_{#pi};counts", tit.Data()));
    /*
    while(hb->GetBinContent(hb->GetMaximumBin())<20){
      hb->RebinX(2);
      if(hb->GetNbinsX()<5)
        break;
    }
    */
    hb->SetMaximum(1.4*hb->GetBinContent(hb->GetMaximumBin()));

    TH1D * hfit=new TH1D(Form("hfit%d", ib),"",1000, 0, 200);
    hfit->SetLineColor(kBlack);
    hfit->SetLineWidth(2);

    Int_t idmax = 0;
    Int_t idmin = 0;

    TH1D *hpar[]={0x0,0x0,0x0,0x0};
    for(Int_t itype=0; itype<4; itype++){
      //pPb only 0-resolution seen
      if(fitMean[itype]<0 || fitReso[itype]<EPSILON){
        printf("continue fit Mean or res negative value! itype %d ib %d %e %e\n", itype, ib, fitMean[itype], fitReso[itype]);
        continue;
      }

      if(fitMean[itype]>fitMean[idmax]){
        idmax = itype;
      }
      if(fitMean[itype]<fitMean[idmin]){
        idmin = itype;
      }

      hpar[itype] = new TH1D(Form("hpar%d%d",ib, itype),"",1000,0,200);
      hpar[itype]->SetLineWidth(2);

      for(Int_t ifit=1; ifit<=hfit->GetNbinsX(); ifit++){
        const Double_t kk = hfit->GetXaxis()->GetBinCenter(ifit);
        const Double_t oldcount = hfit->GetBinContent(ifit);
        const Double_t vfit = hyield->GetBinContent(ib)*fitFrac[itype]*TMath::Gaus(kk, fitMean[itype], fitMean[itype]*fitReso[itype], kTRUE)*hb->GetXaxis()->GetBinWidth(1);
        hpar[itype]->SetBinContent(ifit, vfit);
        hfit->SetBinContent(ifit, oldcount+ vfit);
        //printf("test ib %d itype %d ifit %d kk %e old %e new %e f1 %e frac %e mean %e reso %e\n", ib, itype, ifit, kk, oldcount, hfit->GetBinContent(ifit), vfit, fitFrac[itype], fitMean[itype], fitReso[itype]);
      }

      //printf("test bw %e\n", hb->GetXaxis()->GetBinWidth(ib));
      hpar[itype]->SetLineColor(col[itype]);
    }
    printf("test pp %e idmin %d idmax %d\n", pp, idmin, idmax);

    
    TLegend *lg = new TLegend(0.15,0.5,0.5,0.85);
    style::ResetStyle(lg);
    lg->SetHeader(Form("#splitline{TPC Coherent Fit #it{p}^{track} %.2f GeV/#it{c}}{ d#it{E}/d#it{x}_{#pi} %.1f#pm%.1f (a.u.) #sigma_{#pi} %.1f#pm%.1f%%}", pp, fitMean[1],fitErrMean[1], fitReso[1]*100, fitErrReso[1]*100));
    //hb->SetTitle(Form("%s, #it{p}^{track} %.1f GeV/#it{c};(d#it{E}/d#it{x}-d#it{E}/d#it{x}_{#pi})/#sigma_{#pi};counts", tit.Data(), pp));   

    //hb->GetXaxis()->SetRangeUser(fitMean[idmin]*(1-4*fitReso[idmin]), fitMean[idmax]*(1+4*fitReso[idmax]));
    const Double_t nsig = 15;
    hb->GetXaxis()->SetRangeUser(fitMean[1]*(1-nsig*fitReso[1]), fitMean[1]*(1+nsig*fitReso[1]));
    //hb->GetXaxis()->SetLimits((hb->GetXaxis()->GetXmin()-fitMean[1])/fitMean[1]/fitReso[1], (hb->GetXaxis()->GetXmax()-fitMean[1])/fitMean[1]/fitReso[1]);
    c1->SetRightMargin(0.03);

    //hb->GetXaxis()->SetRangeUser(-10.1,10);
    //hb->GetXaxis()->SetRangeUser(-10.1,15);

    hb->Draw("hist");

//================
    //fit p 0.22 - 0.4
    Double_t pionbias = -999;
    Double_t elebias = -999;
    TF1 *dg = 0x0;
    if(pp>0.16 && pp<0.45){
      dg = new TF1(Form("dg%d",ib),"[0]*(TMath::Gaus(x,[1],[1]*[2],1)+TMath::Abs([3])*TMath::Gaus(x,[1]+TMath::Abs([4]),([1]+TMath::Abs([4]))*[5],1))", 0,200);
      dg->SetParameters(hb->Integral(0,10000),fitMean[1],fitReso[1], fitFrac[3], fitMean[3]-fitMean[1], fitReso[3]);
      hb->Fit(dg,"L");
      Double_t tmppars[10];
      dg->GetParameters(tmppars);
      const Double_t tpion = tmppars[0]; //1./(1+fabs(tmppars[3]));
      const Double_t tele = tmppars[0]*tmppars[3]; //1-tpion;
      const Double_t mpion = fitFrac[1]*hyield->GetBinContent(ib)*hb->GetBinWidth(1);
      const Double_t mele  = fitFrac[3]*hyield->GetBinContent(ib)*hb->GetBinWidth(1);

      pionbias = (mpion-tpion)/tpion*100.;
      elebias = (mele-tele)/tele*100.;
      printf("testing pp %f ib %d fit - true bias pion %f -  %f = %f percent, err %f, electron %f - %f = %f percent, sum fit %f\n", pp, ib, mpion, tpion, pionbias, fitErrFrac[1]/fitFrac[1]*100, mele, tele, elebias, fitFrac[1]+fitFrac[3]);
      dg->SetLineColor(kYellow);
      dg->Draw("same");
    }
    //===============

    const TString parn[]={"p+#bar{p}","#pi^{+}+#pi^{-}","K^{+}+K^{-}","e^{+}+e^{-}"};
    for(Int_t itype=0; itype<4; itype++){
      if(hpar[itype]){
        //hpar[itype]->GetXaxis()->SetLimits((hpar[itype]->GetXaxis()->GetXmin()-fitMean[1])/fitMean[1]/fitReso[1], (hpar[itype]->GetXaxis()->GetXmax()-fitMean[1])/fitMean[1]/fitReso[1]);
        hpar[itype]->Draw("same hist");
        lg->AddEntry(hpar[itype], parn[itype],"l");
      }
    }
    if(dg){
      //lg->AddEntry(dg,Form("#pi bias %.1f%%=%.1f#sigma_{sys.}",pionbias, pionbias/fitErrFrac[1]/100.),"l");
      lg->AddEntry(dg,Form("#pi bias %.1f%%",pionbias),"l");
    }

    hfit->SetLineStyle(kDotted);
    //hfit->GetXaxis()->SetLimits((hfit->GetXaxis()->GetXmin()-fitMean[1])/fitMean[1]/fitReso[1], (hfit->GetXaxis()->GetXmax()-fitMean[1])/fitMean[1]/fitReso[1]);
    hfit->Draw("same hist");
    lg->AddEntry(hfit, "sum","l");

    lg->Draw();

    //gPad->SetGrid(1,0);
    c1->Print(Form("h%02d.eps", ib));
    c1->Print(Form("h%02d.root", ib));

   
  }
}

int main()
{
  /*
see TError.h and set gErrorIgnoreLevel to one of the values:
> const Int_t kPrint = 0;
> const Int_t kInfo = 1000;
> const Int_t kWarning = 2000;
> const Int_t kError = 3000;
> const Int_t kBreak = 4000;
> const Int_t kSysError = 5000;
> const Int_t kFatal = 6000;
   */
  gErrorIgnoreLevel = kWarning;

  style::SetGlobalStyle(0);

  drawMeanGrSys01();
  return 0;
}
